home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 1.iso / HENSA / MATHS / PLPLOT / PLPLOT.ZIP / src / plshade.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-07-13  |  22.5 KB  |  920 lines
  1. /* $Id: plshade.c,v 1.10 1994/07/12 19:20:55 mjl Exp $
  2.  * $Log: plshade.c,v $
  3.  * Revision 1.10  1994/07/12  19:20:55  mjl
  4.  * Argument list for plshade() fixed.
  5.  *
  6.  * Revision 1.9  1994/06/30  18:22:16  mjl
  7.  * All core source files: made another pass to eliminate warnings when using
  8.  * gcc -Wall.  Lots of cleaning up: got rid of includes of math.h or string.h
  9.  * (now included by plplot.h), and other minor changes.  Now each file has
  10.  * global access to the plstream pointer via extern; many accessor functions
  11.  * eliminated as a result.
  12.  *
  13.  * Revision 1.8  1994/04/30  16:15:12  mjl
  14.  * Fixed format field (%ld instead of %d) or introduced casts where
  15.  * appropriate to eliminate warnings given by gcc -Wall.
  16.  *
  17.  * Revision 1.7  1994/03/23  08:29:58  mjl
  18.  * Functions from plctest.c moved to this file for clarity and simplicity.
  19.  * Main shade function now accepts a function evaluator and data much like
  20.  * the contour functions.  Front-ends exist (plshade1(), plshade2()) to call
  21.  * it in a simpler fashion, using predefined memory organizations for the
  22.  * array.  Also, main shade function now accepts two arguments for handling
  23.  * color of the shaded region -- the color map (0 or 1) and the index
  24.  * (a float; >1 and integral for cmap0, in range [0,1] for cmap1).
  25. */
  26.  
  27. /*    plshade.c
  28.  
  29.     Functions to shade regions on the basis of value.
  30.     Can be used to shade contour plots or alone.
  31.     Copyright 1993 Wesley Ebisuzaki 
  32. */
  33.  
  34. /*----------------------------------------------------------------------*\
  35.  * Call syntax for plshade():
  36.  * 
  37.  * void plshade(PLFLT *a, PLINT nx, PLINT ny, char *defined, 
  38.  *    PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax, 
  39.  *     PLFLT shade_min, PLFLT shade_max, 
  40.  *     PLINT sh_color, PLINT sh_width, PLINT min_color, PLINT min_width,
  41.  *     PLINT max_color, PLINT max_width, void (*fill)(), PLINT
  42.  *     rectangular, ...)
  43.  * 
  44.  * arguments:
  45.  * 
  46.  *     PLFLT &(a[0][0])
  47.  * 
  48.  * Contains array to be plotted. The array must have been declared as
  49.  * PLFLT a[nx][ny].  See following note on fortran-style arrays.
  50.  * 
  51.  *     PLINT nx, ny
  52.  * 
  53.  * Dimension of array "a".
  54.  * 
  55.  *     char &(defined[0][0])
  56.  * 
  57.  * Contains array of flags, 1 = data is valid, 0 = data is not valid.
  58.  * This array determines which sections of the data is to be plotted.
  59.  * This argument can be NULL if all the values are valid.  Must have been
  60.  * declared as char defined[nx][ny].
  61.  * 
  62.  *     PLFLT xmin, xmax, ymin, ymax
  63.  * 
  64.  * Defines the "grid" coordinates.  The data a[0][0] has a position of
  65.  * (xmin,ymin).
  66.  * 
  67.  *     void (*mapform)()
  68.  * 
  69.  * Transformation from `grid' coordinates to world coordinates.  This
  70.  * pointer to a function can be NULL in which case the grid coordinates
  71.  * are the same as the world coordinates.
  72.  * 
  73.  *     PLFLT shade_min, shade_max
  74.  * 
  75.  * Defines the interval to be shaded. If shade_max <= shade_min, plshade
  76.  * does nothing.
  77.  * 
  78.  *    PLINT sh_cmap, PLFLT sh_color, PLINT sh_width
  79.  * 
  80.  * Defines color map, color map index, and width used by the fill pattern.
  81.  * 
  82.  *     PLINT min_color, min_width, max_color, max_width
  83.  * 
  84.  * Defines pen color, width used by the boundary of shaded region. The min
  85.  * values are used for the shade_min boundary, and the max values are used
  86.  * on the shade_max boundary.  Set color and width to zero for no plotted
  87.  * boundaries.
  88.  * 
  89.  *     void (*fill)()
  90.  * 
  91.  * Routine used to fill the region.  Use plfill.  Future version of plplot
  92.  * may have other fill routines.
  93.  * 
  94.  *     PLINT rectangular
  95.  * 
  96.  * Flag. Set to 1 if rectangles map to rectangles after (*mapform)() else
  97.  * set to zero. If rectangular is set to 1, plshade tries to save time by
  98.  * filling large rectangles.  This optimization fails if (*mapform)()
  99.  * distorts the shape of rectangles.  For example a plot in polor
  100.  * coordinates has to have rectangular set to zero.
  101.  * 
  102.  * Example mapform's:
  103.  * 
  104.  * Grid to world coordinate transformation.
  105.  * This example goes from a r-theta to x-y for a polar plot.
  106.  *
  107.  * void mapform(PLINT n, PLFLT *x, PLFLT *y) {
  108.  *     int i;
  109.  *     double r, theta;
  110.  *     for (i = 0; i < n; i++) {
  111.  *         r = x[i];
  112.  *         theta = y[i];
  113.  *         x[i] = r*cos(theta);
  114.  *         y[i] = r*sin(theta);    
  115.  *     }
  116.  * }
  117.  * 
  118.  * Grid was in cm, convert to world coordinates in inches.
  119.  * Expands in x direction.
  120.  *
  121.  * void mapform(PLINT n, PLFLT *x, PLFLT *y) {
  122.  *     int i;
  123.  *     for (i = 0; i < n; i++) {
  124.  *         x[i] = (1.0 / 2.5) * x[i];
  125.  *         y[i] = (1.0 / 2.5) * y[i];
  126.  *     }
  127.  * }
  128.  *
  129. \*----------------------------------------------------------------------*/
  130.  
  131. #include "plplotP.h"
  132. #include <float.h>
  133.  
  134. #define NEG  1
  135. #define POS  8
  136. #define OK   0
  137. #define UNDEF 64
  138.  
  139. #define linear(val1, val2, level)  ((level - val1) / (val2 - val1))
  140.  
  141. /* Global variables */
  142.  
  143. static PLFLT sh_max, sh_min;
  144. static int min_points, max_points, n_point;
  145. static int min_pts[4], max_pts[4];
  146. static PLINT pen_col_min, pen_col_max;
  147. static PLINT pen_wd_min, pen_wd_max;
  148.  
  149. /* Function prototypes */
  150.  
  151. static void 
  152. set_cond(register int *cond, register PLFLT *a,
  153.      register const char *defined, register PLINT n);
  154.  
  155. static int 
  156. find_interval(PLFLT a0, PLFLT a1, PLINT c0, PLINT c1, PLFLT *x);
  157.  
  158. static void 
  159. poly(void (*fill) (PLINT, PLFLT *, PLFLT *),
  160.      PLFLT *x, PLFLT *y, PLINT v1, PLINT v2, PLINT v3, PLINT v4);
  161.  
  162. static void 
  163. big_recl(int *cond_code, register int ny, int dx, int dy,
  164.      int *ix, int *iy);
  165.  
  166. static void 
  167. draw_boundary(PLINT slope, PLFLT *x, PLFLT *y);
  168.  
  169. static PLINT 
  170. plctest(PLFLT *x, PLFLT level);
  171.  
  172. static PLINT 
  173. plctestez(PLFLT *a, PLINT nx, PLINT ny, PLINT ix,
  174.       PLINT iy, PLFLT level);
  175.  
  176. /*----------------------------------------------------------------------*\
  177.  * plshade()
  178.  *
  179.  * Shade region.
  180.  * This interface to plfshade() assumes the 2d function array is passed
  181.  * via a (PLFLT **), and is column-dominant (normal C ordering).
  182. \*----------------------------------------------------------------------*/
  183.  
  184. void 
  185. c_plshade(PLFLT **a, PLINT nx, PLINT ny, const char **defined,
  186.       PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
  187.       PLFLT shade_min, PLFLT shade_max,
  188.       PLINT sh_cmap, PLFLT sh_color, PLINT sh_width,
  189.       PLINT min_color, PLINT min_width,
  190.       PLINT max_color, PLINT max_width,
  191.       void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
  192.       void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
  193.       PLPointer pltr_data)
  194. {
  195.     PLfGrid2 grid;
  196.  
  197.     grid.f = a;
  198.     grid.nx = nx;
  199.     grid.ny = ny;
  200.  
  201.     plfshade(plf2eval2, (PLPointer) &grid,
  202.          NULL, NULL,
  203. /*         plc2eval, (PLPointer) &cgrid,*/
  204.          nx, ny, xmin, xmax, ymin, ymax, shade_min, shade_max,
  205.          sh_cmap, sh_color, sh_width,
  206.          min_color, min_width, max_color, max_width,
  207.          fill, rectangular, pltr, pltr_data);
  208. }
  209.  
  210. /*----------------------------------------------------------------------*\
  211.  * plshade1()
  212.  *
  213.  * Shade region.
  214.  * This interface to plfshade() assumes the 2d function array is passed
  215.  * via a (PLFLT *), and is column-dominant (normal C ordering).
  216. \*----------------------------------------------------------------------*/
  217.  
  218. void 
  219. plshade1(PLFLT *a, PLINT nx, PLINT ny, const char *defined,
  220.      PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
  221.      PLFLT shade_min, PLFLT shade_max,
  222.      PLINT sh_cmap, PLFLT sh_color, PLINT sh_width,
  223.      PLINT min_color, PLINT min_width,
  224.      PLINT max_color, PLINT max_width,
  225.      void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
  226.      void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
  227.      PLPointer pltr_data)
  228. {
  229.     PLfGrid grid;
  230.  
  231.     grid.f = a;
  232.     grid.nx = nx;
  233.     grid.ny = ny;
  234.  
  235.     plfshade(plf2eval, (PLPointer) &grid,
  236.          NULL, NULL,
  237. /*         plc2eval, (PLPointer) &cgrid,*/
  238.          nx, ny, xmin, xmax, ymin, ymax, shade_min, shade_max,
  239.          sh_cmap, sh_color, sh_width,
  240.          min_color, min_width, max_color, max_width,
  241.          fill, rectangular, pltr, pltr_data);
  242. }
  243.  
  244. /*----------------------------------------------------------------------*\
  245.  * plfshade()
  246.  *
  247.  * Shade region.  
  248.  * Array values are determined by the input function and the passed data.
  249. \*----------------------------------------------------------------------*/
  250.  
  251. void 
  252. plfshade(PLFLT (*f2eval) (PLINT, PLINT, PLPointer),
  253.      PLPointer f2eval_data,
  254.      PLFLT (*c2eval) (PLINT, PLINT, PLPointer),
  255.      PLPointer c2eval_data,
  256.      PLINT nx, PLINT ny, 
  257.      PLFLT xmin, PLFLT xmax, PLFLT ymin, PLFLT ymax,
  258.      PLFLT shade_min, PLFLT shade_max,
  259.      PLINT sh_cmap, PLFLT sh_color, PLINT sh_width,
  260.      PLINT min_color, PLINT min_width,
  261.      PLINT max_color, PLINT max_width,
  262.      void (*fill) (PLINT, PLFLT *, PLFLT *), PLINT rectangular,
  263.      void (*pltr) (PLFLT, PLFLT, PLFLT *, PLFLT *, PLPointer),
  264.      PLPointer pltr_data)
  265. {
  266.     PLINT init_width, n, slope, ix, iy;
  267.     int count, i, j;
  268.     PLFLT *a, *a0, *a1, dx, dy;
  269.     PLFLT x[8], y[8], xp[2], tx, ty;
  270.  
  271.     int *c, *c0, *c1;
  272.  
  273.     if (plsc->level < 3) {
  274.     plabort("plfshade: window must be set up first");
  275.     return;
  276.     }
  277.  
  278.     if (nx <= 0 || ny <= 0) {
  279.     plabort("plfshade: nx and ny must be positive");
  280.     return;
  281.     }
  282.  
  283.     if (shade_min >= shade_max) {
  284.     plabort("plfshade: shade_max must exceed shade_min");
  285.     return;
  286.     }
  287.  
  288.     if (pltr == NULL)
  289.     rectangular = 1;
  290.  
  291.     init_width = plsc->width;
  292.  
  293.     pen_col_min = min_color;
  294.     pen_col_max = max_color;
  295.  
  296.     pen_wd_min = min_width;
  297.     pen_wd_max = max_width;
  298.  
  299.     plstyl((PLINT) 0, NULL, NULL);
  300.     plwid(sh_width);
  301.  
  302.     switch (sh_cmap) {
  303.     case 0:
  304.     plcol0((PLINT) sh_color);
  305.     break;
  306.     case 1:
  307.     plcol1(sh_color);
  308.     break;
  309.     default:
  310.     plabort("plfshade: invalid color map selection");
  311.     return;
  312.     }
  313.  
  314.     /* alloc space for value array, and initialize */
  315.     /* This is only a temporary kludge */
  316.  
  317.     if ((a = (PLFLT *) malloc(nx * ny * sizeof(PLFLT))) == NULL) {
  318.     plabort("plfshade: unable to allocate memory for value array");
  319.     return;
  320.     }
  321.  
  322.     for (ix = 0; ix < nx; ix++) 
  323.     for (iy = 0; iy < ny; iy++) 
  324.         a[iy + ix*ny] = f2eval(ix, iy, f2eval_data);
  325.  
  326.     /* alloc space for condition codes */
  327.  
  328.     if ((c = (int *) malloc(nx * ny * sizeof(int))) == NULL) {
  329.     plabort("plfshade: unable to allocate memory for condition codes");
  330.     free(a);
  331.     return;
  332.     }
  333.  
  334.     sh_min = shade_min;
  335.     sh_max = shade_max;
  336.  
  337.     /* Ignore defined array for now */
  338.  
  339.     set_cond(c, a, NULL, nx * ny);
  340.  
  341.     dx = (xmax - xmin) / (nx - 1);
  342.     dy = (ymax - ymin) / (ny - 1);
  343.     a0 = a;
  344.     a1 = a + ny;
  345.     c0 = c;
  346.     c1 = c + ny;
  347.  
  348.     for (ix = 0; ix < nx - 1; ix++) {
  349.  
  350.     for (iy = 0; iy < ny - 1; iy++) {
  351.  
  352.         count = c0[iy] + c0[iy + 1] + c1[iy] + c1[iy + 1];
  353.  
  354.         /* No filling needs to be done for these cases */
  355.  
  356.         if (count >= UNDEF)
  357.         continue;
  358.         if (count == 4 * POS)
  359.         continue;
  360.         if (count == 4 * NEG)
  361.         continue;
  362.  
  363.         /* Entire rectangle can be filled */
  364.  
  365.         if (count == 4 * OK) {
  366.         /* find bigest rectangle that fits */
  367.         if (rectangular) {
  368.             big_recl(c0 + iy, ny, nx - ix, ny - iy, &i, &j);
  369.         }
  370.         else {
  371.             i = j = 1;
  372.         }
  373.         x[0] = x[1] = ix;
  374.         x[2] = x[3] = ix+i;
  375.         y[0] = y[3] = iy;
  376.         y[1] = y[2] = iy+j;
  377.  
  378.         if (pltr) {
  379.             for (i = 0; i < 4; i++) {
  380.             (*pltr) (x[i], y[i], &tx, &ty, pltr_data);
  381.             x[i] = tx;
  382.             y[i] = ty;
  383.             }
  384.         }
  385.         else {
  386.             for (i = 0; i < 4; i++) {
  387.             x[i] = xmin + x[i]*dx;
  388.             y[i] = ymin + y[i]*dy;
  389.             }
  390.         }
  391.         if (fill)
  392.             (*fill) ((PLINT) 4, x, y);
  393.         iy += j - 1;
  394.         continue;
  395.         }
  396.  
  397.         /* Only part of rectangle can be filled */
  398.  
  399.         n_point = min_points = max_points = 0;
  400.         n = find_interval(a0[iy], a0[iy + 1], c0[iy], c0[iy + 1], xp);
  401.         for (j = 0; j < n; j++) {
  402.         x[j] = ix;
  403.         y[j] = iy + xp[j];
  404.         }
  405.  
  406.         i = find_interval(a0[iy + 1], a1[iy + 1],
  407.                   c0[iy + 1], c1[iy + 1], xp);
  408.  
  409.         for (j = 0; j < i; j++) {
  410.         x[j + n] = ix + xp[j];
  411.         y[j + n] = iy + 1;
  412.         }
  413.         n += i;
  414.  
  415.         i = find_interval(a1[iy + 1], a1[iy], c1[iy + 1], c1[iy], xp);
  416.         for (j = 0; j < i; j++) {
  417.         x[n + j] = ix + 1;
  418.         y[n + j] = iy + 1 - xp[j];
  419.         }
  420.         n += i;
  421.  
  422.         i = find_interval(a1[iy], a0[iy], c1[iy], c0[iy], xp);
  423.         for (j = 0; j < i; j++) {
  424.         x[n + j] = ix + 1 - xp[j];
  425.         y[n + j] = iy;
  426.         }
  427.         n += i;
  428.  
  429.         if (pltr) {
  430.         for (i = 0; i < n; i++) {
  431.             (*pltr) (x[i], y[i], &tx, &ty, pltr_data);
  432.             x[i] = tx;
  433.             y[i] = ty;
  434.         }
  435.         }
  436.         else {
  437.         for (i = 0; i < n; i++) {
  438.             x[i] = xmin + x[i]*dx;
  439.             y[i] = ymin + y[i]*dy;
  440.         }
  441.         }
  442.  
  443.         if (min_points == 4)
  444.         slope = plctestez(a, nx, ny, ix, iy, shade_min);
  445.         if (max_points == 4)
  446.         slope = plctestez(a, nx, ny, ix, iy, shade_max);
  447.  
  448.         /* special cases: check number of times a contour is in a box */
  449.  
  450.         switch ((min_points << 3) + max_points) {
  451.           case 000:
  452.           case 020:
  453.           case 002:
  454.           case 022:
  455.         if (fill)
  456.             (*fill) (n, x, y);
  457.         break;
  458.           case 040:    /* 2 contour lines in box */
  459.           case 004:
  460.         if (n != 6)
  461.             fprintf(stderr, "plfshade err n=%d !6", (int) n);
  462.  
  463.         if (slope == 1 && c0[iy] == OK) {
  464.             if (fill)
  465.             (*fill) (n, x, y);
  466.         }
  467.         else if (slope == 1) {
  468.             poly(fill, x, y, 0, 1, 2, -1);
  469.             poly(fill, x, y, 3, 4, 5, -1);
  470.         }
  471.         else if (c0[iy + 1] == OK) {
  472.             if (fill)
  473.             (*fill) (n, x, y);
  474.         }
  475.         else {
  476.             poly(fill, x, y, 0, 1, 5, -1);
  477.             poly(fill, x, y, 2, 3, 4, -1);
  478.         }
  479.         break;
  480.           case 044:
  481.         if (n != 8)
  482.             fprintf(stderr, "plfshade err n=%d !8", (int) n);
  483.         if (fill == NULL)
  484.             break;
  485.         if (slope == 1) {
  486.             poly(fill, x, y, 0, 1, 2, 3);
  487.             poly(fill, x, y, 4, 5, 6, 7);
  488.         }
  489.         else {
  490.             poly(fill, x, y, 0, 1, 6, 7);
  491.             poly(fill, x, y, 2, 3, 4, 5);
  492.         }
  493.         break;
  494.           case 024:
  495.           case 042:
  496.         /* 3 contours */
  497.         if (max_points == 4)
  498.             i = NEG;
  499.         else
  500.             i = POS;
  501.  
  502.         if (c0[iy] == i) {
  503.             slope = NEG;
  504.             poly(fill, x, y, 0, 1, 5, 6);
  505.             poly(fill, x, y, 2, 3, 4, -1);
  506.         }
  507.         else if (c0[iy + 1] == i) {
  508.             slope = POS;
  509.             poly(fill, x, y, 0, 1, 2, 3);
  510.             poly(fill, x, y, 4, 5, 6, -1);
  511.         }
  512.         else if (c1[iy + 1] == i) {
  513.             slope = NEG;
  514.             poly(fill, x, y, 0, 1, 6, -1);
  515.             poly(fill, x, y, 2, 3, 4, 5);
  516.         }
  517.         else if (c1[iy] == i) {
  518.             slope = POS;
  519.             poly(fill, x, y, 0, 1, 2, -1);
  520.             poly(fill, x, y, 3, 4, 5, 6);
  521.         }
  522.         break;
  523.           default:
  524.         fprintf(stderr, "prog err switch\n");
  525.         break;
  526.         }
  527.         draw_boundary(slope, x, y);
  528.  
  529.         plwid(sh_width);
  530.         switch (sh_cmap) {
  531.         case 0:
  532.         plcol0((PLINT) sh_color);
  533.         break;
  534.         case 1:
  535.         plcol1(sh_color);
  536.         break;
  537.         }
  538.     }
  539.  
  540.     a0 = a1;
  541.     c0 = c1;
  542.     a1 += ny;
  543.     c1 += ny;
  544.     }
  545.  
  546.     free(c);
  547.     free(a);
  548.     plwid(init_width);
  549. }
  550.  
  551. /*----------------------------------------------------------------------*\
  552.  * set_cond()
  553.  *
  554.  * Fills out condition code array.
  555. \*----------------------------------------------------------------------*/
  556.  
  557. static void 
  558. set_cond(register int *cond, register PLFLT *a,
  559.      register const char *defined, register PLINT n)
  560. {
  561.     if (defined) {
  562.     while (n--) {
  563.         if (*defined++ == 0)
  564.         *cond++ = UNDEF;
  565.         else if (*a < sh_min)
  566.         *cond++ = NEG;
  567.         else if (*a > sh_max)
  568.         *cond++ = POS;
  569.         else
  570.         *cond++ = OK;
  571.         a++;
  572.     }
  573.     }
  574.     else {
  575.     while (n--) {
  576.         if (*a < sh_min)
  577.         *cond++ = NEG;
  578.         else if (*a > sh_max)
  579.         *cond++ = POS;
  580.         else
  581.         *cond++ = OK;
  582.         a++;
  583.     }
  584.     }
  585. }
  586.  
  587. /*----------------------------------------------------------------------*\
  588.  * find_interval()
  589.  *
  590.  * Two points x(0) = a0, (condition code c0) x(1) = a1, (condition code c1)
  591.  * return interval on the line that are shaded
  592.  * 
  593.  * returns 0 : no points to be shaded 1 : x[0] <= x < 1 is the interval 2 :
  594.  * x[0] <= x <= x[1] < 1 interval to be shaded n_point, max_points,
  595.  * min_points are incremented location of min/max_points are stored 
  596. \*----------------------------------------------------------------------*/
  597.  
  598. static int 
  599. find_interval(PLFLT a0, PLFLT a1, PLINT c0, PLINT c1, PLFLT *x)
  600. {
  601.     register int n;
  602.  
  603.     n = 0;
  604.     if (c0 == OK) {
  605.     x[n++] = 0.0;
  606.     n_point++;
  607.     }
  608.     if (c0 == c1)
  609.     return n;
  610.  
  611.     if (c0 == NEG || c1 == POS) {
  612.     if (c0 == NEG) {
  613.         x[n++] = linear(a0, a1, sh_min);
  614.         min_pts[min_points++] = n_point++;
  615.     }
  616.     if (c1 == POS) {
  617.         x[n++] = linear(a0, a1, sh_max);
  618.         max_pts[max_points++] = n_point++;
  619.     }
  620.     }
  621.     if (c0 == POS || c1 == NEG) {
  622.     if (c0 == POS) {
  623.         x[n++] = linear(a0, a1, sh_max);
  624.         max_pts[max_points++] = n_point++;
  625.     }
  626.     if (c1 == NEG) {
  627.         x[n++] = linear(a0, a1, sh_min);
  628.         min_pts[min_points++] = n_point++;
  629.     }
  630.     }
  631.     return n;
  632. }
  633.  
  634. /*----------------------------------------------------------------------*\
  635.  * poly()
  636.  *
  637.  * Draws a polygon from points in x[] and y[].
  638.  * Point selected by v1..v4 
  639. \*----------------------------------------------------------------------*/
  640.  
  641. static void 
  642. poly(void (*fill) (PLINT, PLFLT *, PLFLT *),
  643.      PLFLT *x, PLFLT *y, PLINT v1, PLINT v2, PLINT v3, PLINT v4)
  644. {
  645.     register PLINT n = 0;
  646.     PLFLT xx[4], yy[4];
  647.  
  648.     if (fill == NULL)
  649.     return;
  650.     if (v1 >= 0) {
  651.     xx[n] = x[v1];
  652.     yy[n++] = y[v1];
  653.     }
  654.     if (v2 >= 0) {
  655.     xx[n] = x[v2];
  656.     yy[n++] = y[v2];
  657.     }
  658.     if (v3 >= 0) {
  659.     xx[n] = x[v3];
  660.     yy[n++] = y[v3];
  661.     }
  662.     if (v4 >= 0) {
  663.     xx[n] = x[v4];
  664.     yy[n++] = y[v4];
  665.     }
  666.     (*fill) (n, xx, yy);
  667. }
  668.  
  669. /*----------------------------------------------------------------------*\
  670.  * big_recl()
  671.  *
  672.  * find a big rectangle for shading
  673.  *
  674.  * 2 goals: minimize calls to (*fill)()
  675.  *  keep ratio 1:3 for biggest rectangle
  676.  *
  677.  * only called by plshade()
  678.  *
  679.  * assumed that a 1 x 1 square already fits
  680.  *
  681.  * c[] = condition codes
  682.  * ny = c[1][0] == c[ny]  (you know what I mean)
  683.  *
  684.  * returns ix, iy = length of rectangle in grid units
  685.  *
  686.  * ix < dx - 1
  687.  * iy < dy - 1
  688.  *
  689.  * If iy == 1 -> ix = 1 (so that cond code can be set to skip)
  690. \*----------------------------------------------------------------------*/
  691.  
  692. #define RATIO 3
  693. #define COND(x,y) cond_code[x*ny + y]
  694.  
  695. static void 
  696. big_recl(int *cond_code, register int ny, int dx, int dy,
  697.      int *ix, int *iy)
  698. {
  699.  
  700.     int ok_x, ok_y, j;
  701.     register int i, x, y;
  702.     register int *cond;
  703.  
  704.     /* ok_x = ok to expand in x direction */
  705.     /* x = current number of points in x direction */
  706.  
  707.     ok_x = ok_y = 1;
  708.     x = y = 2;
  709.  
  710.     while (ok_x || ok_y) {
  711. #ifdef RATIO
  712.     if (RATIO * x <= y || RATIO * y <= x)
  713.         break;
  714. #endif
  715.     if (ok_y) {
  716.         /* expand in vertical */
  717.         ok_y = 0;
  718.         if (y == dy)
  719.         continue;
  720.         cond = &COND(0, y);
  721.         for (i = 0; i < x; i++) {
  722.         if (*cond != OK)
  723.             break;
  724.         cond += ny;
  725.         }
  726.         if (i == x) {
  727.         /* row is ok */
  728.         y++;
  729.         ok_y = 1;
  730.         }
  731.     }
  732.     if (ok_x) {
  733.         if (y == 2)
  734.         break;
  735.         /* expand in x direction */
  736.         ok_x = 0;
  737.         if (x == dx)
  738.         continue;
  739.         cond = &COND(x, 0);
  740.         for (i = 0; i < y; i++) {
  741.         if (*cond++ != OK)
  742.             break;
  743.         }
  744.         if (i == y) {
  745.         /* column is OK */
  746.         x++;
  747.         ok_x = 1;
  748.         }
  749.     }
  750.     }
  751.  
  752.     /* found the largest rectangle of 'ix' by 'iy' */
  753.     *ix = --x;
  754.     *iy = --y;
  755.  
  756.     /* set condition code to UNDEF in interior of rectangle */
  757.  
  758.     for (i = 1; i < x; i++) {
  759.     cond = &COND(i, 1);
  760.     for (j = 1; j < y; j++) {
  761.         *cond++ = UNDEF;
  762.     }
  763.     }
  764. }
  765.  
  766. /*----------------------------------------------------------------------*\
  767.  * draw_boundary()
  768.  *
  769.  * Draw boundaries of contour regions based on min_pts[], and max_pts[].
  770. \*----------------------------------------------------------------------*/
  771.  
  772. static void 
  773. draw_boundary(PLINT slope, PLFLT *x, PLFLT *y)
  774. {
  775.     int i;
  776.  
  777.     if (pen_col_min != 0 && pen_wd_min != 0 && min_points != 0) {
  778.     plcol0(pen_col_min);
  779.     plwid(pen_wd_min);
  780.     if (min_points == 4 && slope == 0) {
  781.         /* swap points 1 and 3 */
  782.         i = min_pts[1];
  783.         min_pts[1] = min_pts[3];
  784.         min_pts[3] = i;
  785.     }
  786.     pljoin(x[min_pts[0]], y[min_pts[0]], x[min_pts[1]], y[min_pts[1]]);
  787.     if (min_points == 4) {
  788.         pljoin(x[min_pts[2]], y[min_pts[2]], x[min_pts[3]],
  789.            y[min_pts[3]]);
  790.     }
  791.     }
  792.     if (pen_col_max != 0 && pen_wd_max != 0 && max_points != 0) {
  793.     plcol0(pen_col_max);
  794.     plwid(pen_wd_max);
  795.     if (max_points == 4 && slope == 0) {
  796.         /* swap points 1 and 3 */
  797.         i = max_pts[1];
  798.         max_pts[1] = max_pts[3];
  799.         max_pts[3] = i;
  800.     }
  801.     pljoin(x[max_pts[0]], y[max_pts[0]], x[max_pts[1]], y[max_pts[1]]);
  802.     if (max_points == 4) {
  803.         pljoin(x[max_pts[2]], y[max_pts[2]], x[max_pts[3]],
  804.            y[max_pts[3]]);
  805.     }
  806.     }
  807. }
  808.  
  809. /*----------------------------------------------------------------------*\
  810.  *
  811.  * plctest( &(x[0][0]), PLFLT level)
  812.  * where x was defined as PLFLT x[4][4];
  813.  *
  814.  * determines if the contours associated with level have
  815.  * postive slope or negative slope in the box:
  816.  *
  817.  *  (2,3)   (3,3)
  818.  *
  819.  *  (2,2)   (3,2)
  820.  *
  821.  * this is heuristic and can be changed by the user
  822.  *
  823.  * return 1 if positive slope
  824.  *        0 if negative slope
  825.  *
  826.  * algorithmn:
  827.  *      1st test:
  828.  *    find length of contours assuming positive and negative slopes
  829.  *      if the length of the negative slope contours is much bigger
  830.  *    than the positive slope, then the slope is positive.
  831.  *      (and vice versa)
  832.  *      (this test tries to minimize the length of contours)
  833.  *
  834.  *      2nd test:
  835.  *      if abs((top-right corner) - (bottom left corner)) >
  836.  *       abs((top-left corner) - (bottom right corner)) ) then
  837.  *        return negatiave slope.
  838.  *      (this test tries to keep the slope for different contour levels
  839.  *    the same)
  840. \*----------------------------------------------------------------------*/
  841.  
  842. #define X(a,b) (x[a*4+b])
  843. #define POSITIVE_SLOPE (PLINT) 1
  844. #define NEGATIVE_SLOPE (PLINT) 0
  845. #define RATIO_SQ 6.0
  846.  
  847. static PLINT 
  848. plctest(PLFLT *x, PLFLT level)
  849. {
  850.     double a, b;
  851.     double positive, negative, left, right, top, bottom;
  852.  
  853.     /* find positions of lines */
  854.     /* top = x coor of top intersection */
  855.     /* bottom = x coor of bottom intersection */
  856.     /* left = y coor of left intersection */
  857.     /* right = y coor of right intersection */
  858.  
  859.     left = linear(X(1, 1), X(1, 2), level);
  860.     right = linear(X(2, 1), X(2, 2), level);
  861.     top = linear(X(1, 2), X(2, 2), level);
  862.     bottom = linear(X(1, 1), X(2, 1), level);
  863.  
  864.     /* positive = sq(length of positive contours) */
  865.     /* negative = sq(length of negative contours) */
  866.  
  867.     positive = top * top + (1.0 - left) * (1.0 - left) +
  868.     (1.0 - bottom) * (1.0 - bottom) + right * right;
  869.  
  870.     negative = left * left + bottom * bottom +
  871.     (1.0 - top) * (1.0 - top) + (1.0 - right) * (1.0 - right);
  872. #ifdef DEBUG
  873.     fprintf(stderr, "ctest pos %f neg %f lev %f\n", positive, negative, level);
  874. #endif
  875.     if (RATIO_SQ * positive < negative)
  876.     return POSITIVE_SLOPE;
  877.     if (RATIO_SQ * negative < positive)
  878.     return NEGATIVE_SLOPE;
  879.  
  880.     a = X(1, 2) - X(2, 1);
  881.     b = X(1, 1) - X(2, 2);
  882. #ifdef DEBUG
  883.     fprintf(stderr, "ctest a %f  b %f\n", a, b);
  884. #endif
  885.     if (fabs(a) > fabs(b))
  886.     return NEGATIVE_SLOPE;
  887.     return (PLINT) 0;
  888. }
  889.  
  890. /*----------------------------------------------------------------------*\
  891.  * plctestez
  892.  *
  893.  * second routine - easier to use
  894.  * fills in x[4][4] and calls plctest
  895.  *
  896.  * test location a[ix][iy] (lower left corner)
  897. \*----------------------------------------------------------------------*/
  898.  
  899. static PLINT 
  900. plctestez(PLFLT *a, PLINT nx, PLINT ny, PLINT ix,
  901.       PLINT iy, PLFLT level)
  902. {
  903.  
  904.     PLFLT x[4][4];
  905.     int i, j, ii, jj;
  906.  
  907.     for (i = 0; i < 4; i++) {
  908.     ii = ix + i - 1;
  909.     ii = MAX(0, ii);
  910.     ii = MIN(ii, nx - 1);
  911.     for (j = 0; j < 4; j++) {
  912.         jj = iy + j - 1;
  913.         jj = MAX(0, jj);
  914.         jj = MIN(jj, ny - 1);
  915.         x[i][j] = a[ii * ny + jj];
  916.     }
  917.     }
  918.     return plctest(&(x[0][0]), level);
  919. }
  920.